package org.ainlolcat.ainscience.statistics.handlers.fitting.functions;

import org.apache.commons.math.optimization.fitting.ParametricRealFunction;

/**
* @author ainlolcat
*         Date: 5/9/13
*/
public class PsdVoigtFunction implements ParametricRealFunction {

    //y0 + A [ m_U*(2/Pi)*(w_L / (4*dx^2+w_L^2)) + (1-m_U) * sqrt(4ln2/Pi) * 1/w_G * exp((-4ln2/w_G^2)dx^2)]
    //0 - y0
    //1 - A
    //2 - m
    //3 - w_L
    //4 - w_G
    //5 - x0
    public double value(double x, double[] parameters) {
        double y0   = parameters[0];
        double A    = parameters[1];
        double m_U  = parameters[2];
        double w_L  = parameters[3];
        double w_G  = parameters[4];
        double x_C  = parameters[5];
        double dx   = x - x_C;
        double SQRT_4_LN2_OVER_PI =  Math.sqrt(4 * Math.log(2) / Math.PI);                  //sqrt( (4 ln(2))/ Pi )
        double OVER_4DX2_PLUS_W_L2 = 1 / ( 4 * Math.pow(dx,2) + Math.pow(w_L,2) );          //1 / (4*dx^2+w_L^2)
        double EXP  = Math.exp( - 4 * Math.log(2) * Math.pow(dx,2) / Math.pow(w_G,2) );     //exp((-4ln2/w_G^2)dx^2)
        return  y0 +
                A *
                        (
                                m_U * ( 2 / Math.PI ) * w_L * OVER_4DX2_PLUS_W_L2
                                        +
                                        (1 - m_U) * SQRT_4_LN2_OVER_PI *  ( 1 / w_G ) * EXP
                        );
    }

    public double[] gradient(double x, double[] doubles) {
        double[] gradientVector = new double[doubles.length];
        double y0   = doubles[0];
        double A    = doubles[1];
        double m_U  = doubles[2];
        double w_L  = doubles[3];
        double w_G  = doubles[4];
        double x_C  = doubles[5];

        double dx   = x - x_C;
        double OVER_4DX2_PLUS_W_L2 = 1 / ( 4 * Math.pow(dx,2) + Math.pow(w_L,2) );          //1 / (4*dx^2+w_L^2)
        double SQRT_4_LN2_OVER_PI =  Math.sqrt(4 * Math.log(2) / Math.PI);                  //sqrt( (4 ln(2))/ Pi )
        double EXP  = Math.exp( - 4 * Math.log(2) * Math.pow(dx,2) / Math.pow(w_G,2) );     //exp((-4ln2/w_G^2)dx^2)

        gradientVector[0] = 1;                                                              //dF/dy0
        gradientVector[1] =                                                                 //dF/dA
                m_U * ( 2 / Math.PI ) * w_L * OVER_4DX2_PLUS_W_L2
                        +
                        (1 - m_U) * SQRT_4_LN2_OVER_PI *  ( 1 / w_G ) * EXP;
        gradientVector[2] =                                                                 //dF/dm_U
                A *
                        (
                                ( 2 / Math.PI ) * w_L * OVER_4DX2_PLUS_W_L2
                                        -
                                        SQRT_4_LN2_OVER_PI *  ( 1 / w_G ) * EXP
                        );
        gradientVector[3] =                                                                 //dF/dw_L
                A * m_U * (2/Math.PI) * ( ( 1/OVER_4DX2_PLUS_W_L2 - w_L * 2 * w_L) * Math.pow(OVER_4DX2_PLUS_W_L2,2));
        gradientVector[4] =                                                                 //dF/dw_G
                A * ( 1 - m_U ) * SQRT_4_LN2_OVER_PI * ( -1/Math.pow(w_G,2) + (1/w_G)*(-4)*(-2)*Math.log(2)/Math.pow(w_G,3)) * EXP;
        gradientVector[5] =                                                                 //dF/dx_C
                A *
                        (
                                m_U * (2/Math.PI) * w_L *(- 4 * (-2*dx) * Math.pow(OVER_4DX2_PLUS_W_L2,2))
                                        +
                                        ( 1 - m_U ) * SQRT_4_LN2_OVER_PI * ( 1 / w_G )
                                                *
                                                (- 4 * Math.log(2) / Math.pow(w_G,2)) * ( - 2 * dx ) * EXP
                        );
        return gradientVector;
    }
}
